# from scipy.optimize import linprog
# c=[110, 120, 130, 110, 115,-150]    #目标向量
# A =[[1,1,0,0,0, 0],[0,0,1,1,1,0],[8.8,6.1,2.0,4.2,5.0,-6],[-8.8,-6.1,-2.0,-4.2,-5.0,3]]
# b=[[200],[250],[0],[0]]; Aeq=[[1,1,1,1,1,-1]]; beq=[0]
# res=linprog(c,A,b,Aeq,beq,)
# # print("目标函数的最小值：",res.fun)
# # print("最优解为：",res.x)
# print("求解结果如下", res)



# from scipy.optimize import linprog
# c=[-90, -64]    #目标向量
# A =[[1, 1],[12, 8]]; b=[[50],[480]]
# bound=((0,100/3.0),(0,None))
# res=linprog(c,A,b,None,None,bound,method='simplex',options={"disp": True})
# print("求解结果如下：",res)



# import numpy as np
# from matplotlib.pyplot import plot, show, rc, legend, subplot
# from scipy.optimize import curve_fit
# rc('font', size=15)
# rc('font', family='SimHei');t0=np.arange(1, 7)
# x0 = np.array([5.081, 4.611, 5.1177, 9.3775, 11.0574, 11.0524])
# xt = np.polyfit(t0,x0,1); xh1=np.polyval(xt,t0)  #计算预测值
# delta1 = abs((xh1-x0))/x0   #计算相对误差
# x1 = np.cumsum(x0)
# xh2 = lambda t, a, b, c: a * np.exp(b * t) + c
# para, cov = curve_fit(xh2, t0, x1)
# xh21 = xh2(t0, *para) #计算累加数列的预测值
# xh22 = np.r_[xh21[0], np.diff(xh21)]  #计算预测值
# delta2 = np.abs((xh22-x0)/x0)  #计算相对误差
# print("拟合的参数值为：", para)
# subplot(121)
# plot(t0, x0, 's')
# plot(t0, xh1, '*-')
# legend(('原始数据点', '线性拟合'),loc='upper left')
# subplot(122)
# plot(t0, x1, 'o')
# plot(t0, xh21, 'p-')
# legend(('累加数据点', '累加后拟合'))
# show()


import pandas as pd
import numpy as np
from matplotlib.pyplot import plot,show,rc,legend,xticks

def step_ratio(x0):
    n = len(x0)
    ratio = [x0[i]/x0[i+1] for i in range(len(x0)-1)]
    print(f"级比：{ratio}")
    min_la, max_la = min(ratio), max(ratio)
    thred_la = [np.exp(-2/(n+2)), np.exp(2/(n+2))]
    if min_la < thred_la[0] or max_la > thred_la[-1]:
        print("级比超过灰色模型的范围")
    else:
        print("级比满足要求，可用GM(1,1)模型")
    return ratio, thred_la

def predict(x0):
    n = len(x0)
    x1 = np.cumsum(x0)
    z = np.zeros(n-1)
    for i in range(n-1):
        z[i] = 0.5*(x1[i]+x1[i+1])
    B = [-z, [1]*(n-1)]
    Y = x0[1:]
    u = np.dot(np.linalg.inv(np.dot(B, np.transpose(B))),np.dot(B, Y))
    x1_solve = np.zeros(n)
    x0_solve = np.zeros(n)
    x1_solve[0] = x0_solve[0] = x0[0]
    for i in range(1, n):
        x1_solve[i] = (x0[0]-u[1]/u[0])*np.exp(-u[0]*i)+u[1]/u[0]
    for i in range(1, n):
        x0_solve[i] = x1_solve[i] - x1_solve[i-1]
    return x0_solve, x1_solve, u

def accuracy(x0, x0_solve, ratio, u):
    epsilon = x0 - x0_solve
    delta = abs(epsilon / x0)
    print(f"相对误差：{delta}")
    # Q = np.mean(delta)
    # C = np.std(epsilon)/np.std(x0)

    S1 = np.std(x0)
    S1_new = S1*0.6745
    temp_P = epsilon[abs(epsilon-np.mean(epsilon)) < S1_new]
    P = len(temp_P)/len(x0)
    print(f"预测准确率：{P*100}%")
    ratio_solve = [x0_solve[i]/x0_solve[i+1] for i in range(len(x0_solve)-1)]
    rho = [1-(1-0.5*u[0]/u[1])/(1+0.5*u[0]/u[1])*(ratio[i]/ratio_solve[i]) for i in range(len(ratio))]
    print(f"级比偏差：{rho}")
    # print()
    return epsilon, delta, rho, P


if __name__ == '__main__':
    data=pd.DataFrame(data={"year":[1986,1987, 1988, 1989, 1990, 1991, 1992], "eqL":[25723,30379,34473,38485,40514,42400,48337]})
    x0 = np.array(data.iloc[:,1])
    ratio, thred_la = step_ratio(x0)
    x0_solve, x1_solve, u = predict(x0)



    epsilon, delta, rho, P = accuracy(x0, x0_solve, ratio, u)


# -*- coding: utf-8 -*-
# @Time : 2022/3/18 14:18
# @Author : Orange
# @File : g_pred.py.py

# from decimal import *
#
#
# class GM11():
#     def __init__(self):
#         self.f = None
#
#     def isUsable(self, X0):
#         '''判断是否通过光滑检验'''
#         X1 = X0.cumsum()
#         rho = [X0[i] / X1[i - 1] for i in range(1, len(X0))]
#         rho_ratio = [rho[i + 1] / rho[i] for i in range(len(rho) - 1)]
#         print("rho:", rho)
#         print("rho_ratio:", rho_ratio)
#         flag = True
#         for i in range(2, len(rho) - 1):
#             if rho[i] > 0.5 or rho[i + 1] / rho[i] >= 1:
#                 flag = False
#         if rho[-1] > 0.5:
#             flag = False
#         if flag:
#             print("数据通过光滑校验")
#         else:
#             print("该数据未通过光滑校验")
#
#         '''判断是否通过级比检验'''
#         lambds = [X0[i - 1] / X0[i] for i in range(1, len(X0))]
#         X_min = np.e ** (-2 / (len(X0) + 1))
#         X_max = np.e ** (2 / (len(X0) + 1))
#         for lambd in lambds:
#             if lambd < X_min or lambd > X_max:
#                 print('该数据未通过级比检验')
#                 return
#         print('该数据通过级比检验')
#
#     def train(self, X0):
#         X1 = X0.cumsum()
#         Z = (np.array([-0.5 * (X1[k - 1] + X1[k]) for k in range(1, len(X1))])).reshape(len(X1) - 1, 1)
#         # 数据矩阵A、B
#         A = (X0[1:]).reshape(len(Z), 1)
#         B = np.hstack((Z, np.ones(len(Z)).reshape(len(Z), 1)))
#         # 求灰参数
#         a, u = np.linalg.inv(np.matmul(B.T, B)).dot(B.T).dot(A)
#         u = Decimal(u[0])
#         a = Decimal(a[0])
#         print("灰参数a：", a, "，灰参数u：", u)
#         self.f = lambda k: (Decimal(X0[0]) - u / a) * np.exp(-a * k) + u / a
#
#     def predict(self, k):
#         X1_hat = [float(self.f(k)) for k in range(k)]
#         X0_hat = np.diff(X1_hat)
#         X0_hat = np.hstack((X1_hat[0], X0_hat))
#         return X0_hat
#
#     def evaluate(self, X0_hat, X0):
#         '''
#         根据后验差比及小误差概率判断预测结果
#         :param X0_hat: 预测结果
#         :return:
#         '''
#         S1 = np.std(X0, ddof=1)  # 原始数据样本标准差
#         S2 = np.std(X0 - X0_hat, ddof=1)  # 残差数据样本标准差
#         C = S2 / S1  # 后验差比
#         Pe = np.mean(X0 - X0_hat)
#         temp = np.abs((X0 - X0_hat - Pe)) < 0.6745 * S1
#         p = np.count_nonzero(temp) / len(X0)  # 计算小误差概率
#         print("原数据样本标准差：", S1)
#         print("残差样本标准差：", S2)
#         print("后验差比：", C)
#         print("小误差概率p：", p)
#
#
# if __name__ == '__main__':
#     import matplotlib.pyplot as plt
#     import numpy as np
#
#     plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
#     plt.rcParams['axes.unicode_minus'] = False  # 步骤二（解决坐标轴负数的负号显示问题）
#
#     # 原始数据X
#     X = np.array(
#         [21.2, 22.7, 24.36, 26.22, 28.18, 30.16, 32.34, 34.72, 37.3, 40.34, 44.08, 47.92, 51.96, 56.02, 60.14,
#          64.58,
#          68.92, 73.36, 78.98, 86.6])
#     # 训练集
#     X_train = X[:int(len(X) * 0.7)]
#     # 测试集
#     X_test = X[int(len(X) * 0.7):]
#
#     model = GM11()
#     model.isUsable(X_train)  # 判断模型可行性
#     model.train(X_train)  # 训练
#     Y_pred = model.predict(len(X))  # 预测
#     Y_train_pred = Y_pred[:len(X_train)]
#     Y_test_pred = Y_pred[len(X_train):]
#     score_test = model.evaluate(Y_test_pred, X_test)  # 评估
#
#     # 可视化
#     plt.grid()
#     plt.plot(np.arange(len(X_train)), X_train, '->')
#     plt.plot(np.arange(len(X_train)), Y_train_pred, '-o')
#     plt.legend(['负荷实际值', '灰色预测模型预测值'])
#     plt.title('训练集')
#     plt.show()
#
#     plt.grid()
#     plt.plot(np.arange(len(X_test)), X_test, '->')
#     plt.plot(np.arange(len(X_test)), Y_test_pred, '-o')
#     plt.legend(['负荷实际值', '灰色预测模型预测值'])
#     plt.title('测试集')
#     plt.show()









